Feedback predictive control approach for processes with time delay in the manipulated variable

ABSTRACT

This invention relates to a feedback predictive controller, systems comprising and methods employing the same. Preferably the feedback predictive controller and/or systems comprising the feedback predictive controller are part of an automatic insulin delivery system. The methods described herein can be used to control blood glucose concentration in a patient with diabetes. Preferably, the insulin delivery system is an artificial pancreas.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims priority under 35 U.S.C. § 119 to Provisional Application U.S. Ser. No. 62/697,798, filed on Jul. 13, 2018 and Provisional Application U.S. Ser. No. 62/701,153 filed on Jul. 20, 2018, both of which are herein incorporated by reference in their entirety including without limitation, the specification, claims, and abstract, as well as any figures, tables, or examples thereof.

TECHNICAL FIELD

This invention relates generally to a control system that utilizes feedback predictive control. More specifically, but not exclusively, the invention relates to an automatic insulin delivery system that employs feedback predictive control.

BACKGROUND OF THE INVENTION

In a person without diabetes, several systems including, for example, the metabolic, endocrine, and cardiovascular systems, function collectively to maintain homeostasis. However, in a person with diabetes, the inherent glucose regulation mechanism is dysfunctional. Glucose levels are affected by the state of the metabolic-physiological-endocrine system and disturbances to that system. Disturbances such as stress, physical activity, hormonal changes, periods of growth, medications, illness/infection and fatigue, that are unknown in many cases, can cause large changes in blood glucose concentration (BGC) for individuals with type 1 diabetes. Their combined effects on BGC are highly complex, highly correlated, and difficult to measure and model accurately. In a normal person, several systems, such as the metabolic, endocrine, and cardiovascular systems function collectively to maintain homeostasis. However, in a person with diabetes their inherent glucose regulation mechanism becomes dysfunctional. If the amount of carbohydrates in a meal was the only variable to affect BGC, management would be much simpler. However, not only are there other disturbances in food, such as fats and proteins, but also in the levels and types of activities as well as physical and emotional stress, among many others. Furthermore, not only does BGC change due to several disturbances, their impact on glucose level is highly correlated, dynamic, and nonlinear, making it difficult to distinguish the effect each input has on BGC.

Insulin therapy involves multiple daily doses of insulin before meals to correct high blood glucose, with the amount either pre-recommended by a physician or decided by the patient on the basis measured BGC and the number of carbohydrates they estimate will be ingested at the time of the meal. This protocol is inconvenient and unreliable. It often results in hypoglycemic and hyperglycemic episodes, both of which can be life-limiting and life-threatening. Thus, there has been a desire to develop more optimized insulin delivery systems and protocols. This has included research to develop automatic insulin delivery systems.

The potential for successful automatic insulin delivery has entered a new era due to recent technological advancements of insulin pumps and blood glucose sensors. Consequently, for full automation and control capable of reducing the variance in BGC, the control algorithm must be capable of tight control for major disturbances such as meals, activity, and stress. However, the control systems in place to date are not sufficient to address the complexities of BGC control.

For example, classical time-delayed feedback control (FBC) approaches, such as the Smith Predictor, address time delay (θ) in the controlled variable (CV) when a change in the manipulated variable (MV) affects CV immediately but the measurement of CV corresponding to the MV change is much later (i.e., delayed). In contrast, when MV has time delay, a change in MV does not affect CV until θ time later. Thus, in this case, optimal feedback control requires knowledge of the value(s) of CV after it begins to change from the current change in MV, i.e., future predictions of CV a time distant into the future at least as great as θ. Hence, an optimal FBC-type of algorithm for addressing dead time in MV requires future prediction of CV to change MV. The classical FBC-type of algorithm for dead time in MV is model predictive control (MPC). The MPC control law is for the future predicted value of CV to be at the set point an integer value of sampling time steps (J) beyond the amount of dead time. Hence, the effectiveness of MPC depends strongly on the accuracy of the future prediction of CV, which decreases as J increases. Moreover, the larger the time lag for MV, the longer it takes to reach the set point and the greater J will be. Thus, if time lag for MV is significant, the more difficult it will be to predict CV and the greater the potential error there will be in that prediction. For this reason, all types of classical time-delay feedback control algorithms are unable to address the complexities of BGC.

Theoretically, the superiority of feedforward control (FFC) over all other control systems is that corrective action can be taken to cancel the effects of disturbances on the control variable (i.e., BGC) proactively. Currently, FFC is not capable of approaching this theoretical limit, i.e., very tight control of BGC, because insulin, the MV, can only decrease BGC. Before it is possible to approach this idealized control, a new hormone needs to be developed and approved by the FDA that can increase BGC. In addition, FFC has not been capable of addressing the complexity of various measured disturbances impacting the regulation of blood glucose levels associated with diabetes due to the limits of sensor technology as was as effective model development of the complex causative relationships of disturbances and the controlled variable. These insufficiencies will have to be overcome before FFC is able to significantly improve control of BGC. Thus, feedforward control, while being theoretically superior, is at this time, incapable of serving as an approved automatic insulin delivery system.

Consequently, there is a need for a system providing automatic delivery of insulin with minimal variability around the desired BGC target using a predictive feedback control approach. More specifically, currently, the most practical challenge is providing an effective FBC-type of predictive control system in the presence of unmeasured disturbances, with a minimal prediction horizon (i.e., prediction a distance B into the future) for the greatest possible prediction accuracy of CV, which leads to the most effective control system possible for BGC.

Accordingly, it is an objective of the claimed invention to provide an methodology that more accurately controls blood glucose levels for FBC-type predictive control systems.

Other objects, advantages and features of the present invention will become apparent from the following specification taken in conjunction with the accompanying figures.

While multiple embodiments are disclosed, still other embodiments of the present invention will become apparent to those skilled in the art from the following detailed description, which shows and describes illustrative embodiments of the invention. Accordingly, the figures and detailed description are to be regarded as illustrative in nature and not restrictive.

BRIEF DESCRIPTION OF THE FIGURES

FIG. 1A is an exemplary block diagram for a feedback predictive controller (FBPC) employing a traditional controller such as a PID controller.

FIG. 1B is an exemplary block diagram for a feedback predictive controller (FBPC) where employing a model predictive controller.

FIG. 2A shows a diagram of a system or method for controlling BGC according to the invention.

FIG. 2B shows a diagram of a system or method for controlling BGC according to the invention.

FIG. 3A shows a block diagram of a traditional Wiener model.

FIG. 3B shows a block diagram of a Wiener/Semi-Coupled model.

FIG. 3C shows a block diagram of a Wiener/Semi-Coupled model with thirteen specific inputs for disturbances on the system.

FIG. 4 shows the model CSTR used in Example 1.

FIG. 5 shows the training input sequences for M, T_(ci) and C_(Ai) from Example 1.

FIG. 6 shows the fit for the training data with Model 2 as given by Eq. 27.

FIG. 7 shows the fit for the training data with Model 1, the input-only model.

FIG. 8A shows the response of the measured reactor temp, T_(m), the variable to be controlled, when the input signal to the valve, M is fixed at its initial value.

FIG. 8B shows the input signal to the valve, M, fixed at its initial value and with changed input variables as if for a control run.

FIGS. 9A and 9B show the T_(m) and M results, respectively, for the FBPC in Example 1.

FIGS. 9C and 9D show the T_(m) and M results, respectively, for the MPC where J=10 in Example 1.

FIGS. 9E and 9F show the T_(m) and M results, respectively, for the MPC where J=15 in Example 1.

FIG. 10 shows a comparison of the FBPC and MPC for BGC control n steps into the future from Example 3.

FIGS. 11A-11E show the graphical control run results for both the MPC and FBPC for the first five (5) test subjects.

FIGS. 12A-12E show the bolus profiles of the insulin feed rates for the first five (5) test subjects both under an MPC and FBPC.

FIG. 12F shows the insulin feed rate of one meal in detail for test subject 5.

FIG. 13 shows an exemplary FBPC scheme as described herein and based on the glucose monitoring/control system exemplified in FIGS. 2A and 2B.

Various embodiments of the present invention are described in detail with reference to the figures. The reference to various embodiments does not limit the scope of the invention. Figures represented herein are not limitations to the various embodiments according to the invention and are presented for exemplary illustration of the invention.

DETAILED DESCRIPTION

The embodiments of this invention are not limited to particular insulin pumps, apparatuses for monitoring insulin and/or blood glucose levels, activity inputs, and stress measuring devices. Moreover, while the Feedback Predictive Controller and systems employing the same are primarily discussed in the context of diabetes and BGC control, the controller is not limited to diabetes and BGC control. Rather it is suitable for other systems, particularly, those having significant lag time. It is further to be understood that all terminology used herein is for the purpose of describing particular embodiments only and is not intended to be limiting in any manner or scope. For example, all units, prefixes, and symbols may be denoted in its SI accepted form.

Numeric ranges recited within the specification are inclusive of the numbers defining the range and include each integer, fraction, and decimal within the defined range. Throughout this disclosure, various aspects of this invention are presented in a range format. It should be understood that the description in range format is merely for convenience and brevity and should not be construed as an inflexible limitation on the scope of the invention. Accordingly, the description of a range should be considered to have specifically disclosed all the possible sub-ranges as well as individual numerical values within that range. For example, description of a range such as from 1 to 6 should be considered to have specifically disclosed sub-ranges such as from 1 to 3, from 1 to 4, from 1 to 5, from 2 to 4, from 2 to 6, from 3 to 6 etc., as well as individual numbers within that range such as 1, 2, 3, 4, 5, and 6, fractions such as 1¾, 3¼, and 4%, and decimals such as 2.3, 4.1, and 5.7. This applies regardless of the breadth of the range.

So that the present invention may be more readily understood, certain terms are first defined. Unless defined otherwise, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which embodiments of the invention pertain. Many methods and materials similar, modified, or equivalent to those described herein can be used in the practice of the embodiments of the present invention without undue experimentation, the preferred materials and methods are described herein. In describing and claiming the embodiments of the present invention, the following terminology will be used in accordance with the definitions set out below.

Definitions

As used herein the singular forms “a,” “an” and “the” can include plural referents unless the content clearly indicates otherwise.

As used herein the term “BGC” refers to blood glucose concentration. It is not limited to any particular measurement unit.

As used herein the term “FBC” refers to feedback control or any synonyms thereof, including for example, “feed-back control.”

As used herein the term “FBE” refers to feedback error or any synonyms thereof and is the deviation of the controlled variable from its set point value.

As used herein the term “FBPC” refers to a feedback predictive control or feedback predictive controller as described in greater detail herein.

As used herein the term “MPC” refers to model predictive control.

As used herein the term “dead time” refers to the amount of time it takes for the controlled variable to start to respond after a change in the manipulated variable.

As used herein the term “lag” refers to a variable not reaching its final value immediately when the variable begins to change due a change in another variable causing it to change. A variable without lag reaches its final value immediately.

As used herein the term “prediction horizon” refers to the distance into the future that the value of the controlled variable must be estimated.

As used herein the term “controlled variable” refers to the output the control system, the purpose of the control system is to keep the controlled variable at the set point.

As used herein the term “manipulated variable” refers to the input variable that is adjusted to keep the controlled variable in the system at a set point.

As used herein the term “PID controller” refers a proportional-integral-derivative controller used in control systems that require continuously modulated control.

Feedback Predictive Controller

The systems and methods described herein comprise a feedback predictive controller (“FBPC”). Preferably, the FBPC is a model-based feedback predictive control system that uses a predictive feedback error (deviation from set point) to provide changes to the system in an effort to nullify the deviation. Surprisingly it has been found that by employing a feedback controller with a predicted FBE, one can reduce, or preferably minimize, deviation of the controlled variable from its set point in cases where manipulated variables have large dead time and time lag. In an aspect of the invention, the set point of a controlled variable is set by the user of the control system. Thus, the set point is determined by the user of the control system.

In the context of BGC and diabetes, accurate modeling of BGC 30 to 60 minutes into the future is critical for model-based feedback predictive control strategies that use the predictive feedback error (deviation from set point) since it takes 30 to 60 minutes for current changes in insulin flow rate to change BGC. To achieve accurate modeling (i.e., estimation) of the future value of the control variable is essential. Preferably, the models will not drift over time, have minimal bias, and provide at least 70%, 75%, 80%, 85%, 90%, 91%, 92%, 93%, 94%, 95%, 96%, 97%, 98%, or 99% accuracy and fit over at least 24 hours, 36 hours, 48 hours, 3 days, 4, days, 5 days, 6 days, 7 days, 2 weeks, 3 weeks, or one month.

The FBPC system preferably comprises a two-stage method that first applies the Wiener/Semi-Coupled method and obtains the input-only Model 1. It then uses the pre-whitening method to model the serial correlation structure of the Model 1 residuals to address unmeasured disturbances and model bias. This is Model 2 and is the predictive model for future BGC using current and past BGC measurements. According to this system, BGC is predicted a distance into the future equal to the given value of dead time for exogenous insulin flow rate. The future feedback error is determined using this value and a conventional PID controller is used to manipulate insulin flowrate in the present. Preferably, FBPC control uses a PID control algorithm, with the minimum prediction horizon, no specified control horizon, and an inherent cause-and-effect relationship between MV and CV that is “built-in” to the PID structure.

The FBPC and systems and methods employing the same have provide improvements and advantages over existing controllers, systems and methods for controlling BGC. For example, one such advantage is that it does not require a cause-and-effect model for the relationship of insulin flow rate on BGC. This is significant as a cause-and-effect model is very difficult to accurately obtain in real-life. Model-based predictive control algorithms such as FFC and MPC require accurate cause-and-effect modeled relationship of insulin flow rate on BGC, a critical limitation of such approaches.

According to the invention as described herein FBPC is similar to MPC in certain respects; however, MPC uses an FBE that is greater than a distance θ into the future (i.e., a greater prediction horizon). Therefore, in applications where the accuracy of CV declines as the prediction horizon increases, which is common, MPC could perform much worse than FBPC when its prediction horizon is significantly greater than θ. In addition, there are other fundamental differences in these two control algorithms that can result in significantly different control performance. MPC uses a fitted mathematical model for the manipulated variable and FBPC uses a form of the PID control algorithm. Thus, MPC performance depends on the accuracy of this model, which can be difficult to obtain in many applications. Furthermore, MPC has one discrete tuning parameter. Therefore, the one-dimensional tuning space of MPC has the advantage of simplicity over the possible three-dimensional one of FBPC. This is also a disadvantage in many applications. FBPC has a tuning space as large as three-dimensions and thus, more geometric directions for finding optimal tuning parameters. In addition, its tuning parameters are continuous and can be set to any value in a range, whereas the one for MPC can only be set at discrete values, in increments that depend on the sampling rate. Furthermore, the prediction horizon for MPC depends on the dynamic lag of MV. The larger this lag, the greater the prediction horizon, and hence the worse prediction accuracy can be. Since the prediction horizon for FBPC is fixed at θ distance in the future, process lag does not affect its prediction horizon.

In a general class, MPC can be any type of control system that uses a model for the control horizon and a model for the prediction horizon. However, as described above, the model for the control horizon must have an accurate cause-and-effect relationship between the inputs, especially MV, and the output, CV. This is usually very difficult to achieve in practice and required by all MPC type of control algorithms. Hence, a critical advantage of FBPC is that it does not use a control horizon and thus, does not need, or use, a model horizon to make changes in MV.

An exemplary block diagram for FBPC is shown in FIG. 1A. The controller G_(c) is preferably a PID controller. FIG. 1B provides an exemplary block diagram mimicking FIG. 1A, but employing a model predictive controller. In both FIGS. 1A and 1B, the output from G_(v) represents adjustments of the MV. In an automatic insulin delivery system, the output from G_(v) would represent varying the insulin flow rate or quantity of insulin administered to a patient.

According to an aspect of the invention, Model 1 can apply a Wiener model or a Wiener/Semi-Coupled model. A traditional Wiener model is shown in FIG. 3A. Preferred Wiener/Semi-Coupled models include those shown in FIGS. 3B and 3C. In the Wiener model (FIG. 3A), the inputs are represented by x_(i) where i=1, . . . , p, which are the measured noninvasive variables or disturbances. Each input has its own linear dynamic block, G_(i), and each dynamic block has an intermediate, unobservable output, v_(i), which represents the independent dynamic response of its corresponding input. All the intermediate v_(i)'s are collected and passed through an unrestricted static gain block, ƒ(V), to produce the final output, y, which represents BGC. If employed in BGC control system, x_(i) represents things like food, activity, and stress, i.e., disturbances that affect BGC and the output variable y is BGC.

As shown in FIGS. 3B and 3C, all of the inputs, x_(i), pass through a linear dynamic block to produce the unobservable dynamic output variables v_(i). Note that, i=A₁, . . . , A_(p), for p activity inputs, C for carbohydrates, F for fats, P for proteins, and I for blood insulin concentration (BIC). FIG. 3C shows an exemplary semi-coupled network 14 accounting for thirteen inputs 16 (disturbances). This includes two insulin inputs, three consumed energy inputs, and eight activity inputs.

In FIGS. 3A-3C, “Y” represents BGC; “V^(set)” denotes the set point for BGC, i.e., the desired value at θ; “x_(a1)” to “x_(I2)” represent measured input variables (e.g. physical activities, insulin, and meal logs); all “v's” are dynamic effects of corresponding “x” inputs; “G_(BI)” represents the BGC due to insulin only; “G_(FI)” represents the BGC due to food and insulin only; “G₁” to “G_(p)” represent dynamic block for each input; “V_(FI)” provides feedback to G_(FI) block of the BGC due to food and insulin changes only.

Consumed Energy

There are three main categories of consumed energy: carbohydrates, fats, and proteins. As shown in FIGS. 3A and 3B, these are accounted for by X_(C), X_(F), and X_(P). The consumed energy outputs feed into the G_(FI) block along with the output from the G_(BI) block. The output from the G_(FI) block then combines with the activity disturbances.

Insulin

There are two types of insulin: basal and bolus. Basal insulin, often referred to as background insulin, functions to keep blood glucose levels consistent during periods of fasting. Basal insulin is produced by the pancreas. Basal insulin can also be provided by the continuous administration of fast-acting insulin in small doses. Bolus insulin is faster acting insulin provided in some dosage. It can be administered intravenously, by intramuscular injection, by intrathecal injection, or by subcutaneous injection. In FIGS. 3B-3C, these are accounted for by X_(I), in particular, X_(I1) and X_(I2). As shown in FIGS. 3B-3C, the insulin inputs combine to form the G_(BI) block, which feeds into the G_(FI) block. Preferably, the G_(FI) block also contains the output from the consumed energy inputs.

Activities

There are many activities that affect BGC. These activities include, but are not limited to, exercise, walking, working, and casual movement. These activities can be monitored by monitoring and/or measuring body position, movement, heat dissipated, skin temperature, near body temperature, galvanic skin response, and sleep, etc. In order to monitor body position and movement, it is preferable to use a two-axis accelerometer. Sleep can be monitored by the lack of certain kinds of movement. The heat dissipated from the body can be accounted for based on the body's heat flux. Skin temperature and near-body temperature can be measured by sensitive thermistors. Galvanic skin response (GSR) can be measured by the conductivity of a person's skin as it varies due to physical and emotional stimuli.

Any type and number of monitoring systems can be used to monitor the various activities. In a preferred embodiment, a single monitoring system is employed to monitor all activities at once. A preferred monitoring system is a type similar to the SenseWear® Pro3 Body Monitoring System, previously available from BodyMedia, Inc. in Pittsburgh, Pa.

Feedback Predictive Control Model

The dynamic blocks for the inputs have second-order-plus-dead-time-plus-lead (SOPDTPL) form as shown in Eq. 1 below:

$\begin{matrix} {{{\tau_{i}^{2}\frac{d^{2}{v_{i}(t)}}{{dt}^{2}}} + {2\tau_{i}\zeta_{i}\frac{{dv}_{i}(t)}{dt}} + {v_{i}(t)}} = {{\tau_{ai}\frac{{dx}_{i}\left( {t - \theta} \right)}{dt}} + {x_{i}\left( {t - \theta_{i}} \right)}}} & (1) \end{matrix}$

where τ_(i) is the time constant, ζ_(i) is the damping coefficient, τ_(ai) is the lead parameter and θ_(i) is the dead time. Using a backward difference approximation

$\left( {{e.g.},{\frac{{dv}_{1}(t)}{dt} \approx \frac{v_{1,t} - v_{1,{t - {\Delta \; t}}}}{\Delta \; t}}} \right)$

applied to a sampling interval of Δt, an approximate discrete-time form of Eq. 1 is:

v _(i,t)=δ_(1,i) v _(i,t−Δt)+δ_(2,i) v _(i,t−2Δt)+ω_(1,i) x _(i,t−k) _(i) _(Δt)+ω_(2,i) x _(i,t−(k) _(i) _(+1)Δt)  (2)

where ω_(2,i)=1−δ_(1,i)−δ_(2,i)−ω_(1,i) to satisfy the unity gain constraint with

$\begin{matrix} {\delta_{1,i} = \frac{{2\tau_{i}^{2}} + {2\tau_{i}\zeta_{i}\Delta \; t}}{\tau_{i}^{2} + {2\tau_{i}\zeta_{i}\Delta \; t} + {\Delta \; t^{2}}}} & (3) \\ {\delta_{2,i} = \frac{- \tau_{i}^{2}}{\tau_{i}^{2} + {2\tau_{i}\zeta_{i}\Delta \; t} + {\Delta \; t^{2}}}} & (4) \\ {\omega_{1,i} = \frac{\left( {\tau_{ai} + {\Delta \; t}} \right)\Delta \; t}{\tau_{i}^{2} + {2\tau_{i}\zeta_{i}\Delta \; t} + {\Delta \; t^{2}}}} & (5) \end{matrix}$

A dynamic mass balance on G_(FI) block in FIGS. 3A-3C, which represents the BGC due to food and insulin only, gives

$\begin{matrix} {\frac{{dV}_{FI}(t)}{dt} = {{a_{C}{v_{C}(t)}} + {a_{F}{v_{F}(t)}} + {a_{P}{v_{P}(t)}} - {a_{FI}{I(t)}{V_{FI}(t)}}}} & (6) \end{matrix}$

where a_(C), a_(F), a_(P), and a_(FI), are estimated model parameters; I(t) is the unmeasured blood-insulin concentration (BIC) at time t; and V_(FI)(t) is the BGC due to food and insulin changes only. Similarly, a dynamic mass balance on the GBI block, which represents BIC, gives

$\begin{matrix} {\frac{{dI}(t)}{dt} = {{a_{I}{v_{I}(t)}} + {a_{Bi}{I(t)}}}} & (7) \end{matrix}$

where a_(I) and a_(BI), are estimated model parameters. The function ƒ(V) is called “the static function.” This function can theoretically be of any form. For effectiveness under mild extrapolation, in modeling real BGC data we will use a first-order linear regression structure, given in discrete form as:

ƒ(V)=η_(t)=ε_(i) =a ₀ +V _(FI,t) +a _(A1) v _(A) _(t) _(,t) + . . . +a _(Ap) v _(A) _(p) _(,t)

where a₀ and a_(Ai)'s are static estimated model parameters.

Preferably, the FBPC controller models unmeasured disturbances and bias in

Model 2. Specifically, a noise model structure is applied to Eq. 8 as follows:

$\begin{matrix} {\mspace{79mu} {{y_{t} = {\eta_{t} + N_{t}}},\mspace{79mu} {where}}} & (9) \\ {N_{t} = {{\frac{\theta_{q}(B)}{\varphi_{p}(B)}a_{t}} = {{\frac{1 - {\theta_{1}B} - {\theta_{2}B^{2}} - \ldots - {\theta_{q}B^{q}}}{1 - {\varphi_{1}B} - {\varphi_{2}B^{2}} - \ldots - {\varphi_{p}B^{p}}}a_{t}} = {\frac{a_{t}}{\Phi (B)} = {\frac{a_{t}}{\left\lbrack {1 - {\phi_{1}B} - {\phi_{2}B^{2}} - \ldots}\mspace{11mu} \right\rbrack}{\forall t}}}}}} & \left( {10a} \right) \\ {\mspace{79mu} {{{subject}\mspace{14mu} {to}\mspace{14mu} {\sum\phi_{i}}} = 1}} & \left( {10b} \right) \\ {\mspace{79mu} {a_{t}\overset{indep}{\sim}{N\left( {0,\sigma^{2}} \right)}}} & (11) \end{matrix}$

such that y_(t) is the measured BGC at t and and B^(γ)x_(t)=X_(t-γΔt). Eq. 10b is a critical innovation of this work as it eliminates model bias, i.e., model drift of the CV over time. Then from substituting Eq. 10 into Eq. 9 and rearranging gives:

γ_(t)└1−φ₁ B−φ ₂ B ²− . . . ┘=η_(t)└1−φ₁ B−φ ₂ B ² − . . . ┘+a _(t)  (12)

With dead time θ=kΔt, Eq. 12 is

$\begin{matrix} {{y_{t}\left\lfloor {1 - {\phi_{1}B^{k}} - {\phi_{2}B^{k + 1}} - \ldots}\mspace{11mu} \right\rfloor} = \left. {{\eta_{t}\left\lfloor {1 - {\phi_{1}B^{k}} - {\phi_{2}B^{k + 1}} - \ldots}\mspace{11mu} \right\rfloor} + a_{t}}\Rightarrow \right.} & (13) \\ {y_{t} = {\eta_{t} + {\phi_{1}\left( {y_{t - {k\; \Delta \; t}} - \eta_{t - {k\; \Delta \; t}}} \right)} + {\phi_{2}\left( {y_{t - {{({k + 1})}\Delta \; t}} - \eta_{t - {{({k + 1})}\Delta \; t}}} \right)} + \ldots + a_{t}}} & (14) \end{matrix}$

Modification of Eq. 14 for predicting θ=kΔt time steps into the future gives

$\begin{matrix} \begin{matrix} {{\hat{y}}_{t + {k\; \Delta \; t}} = {{\hat{\eta}}_{t + {k\; \Delta \; t}} + {{\hat{\phi}}_{1}\left( {y_{t} - {\hat{\eta}}_{t}} \right)} + {{\hat{\phi}}_{2}\left( {{\hat{y}}_{t - {\Delta \; t}} - {\hat{\eta}}_{t - {\Delta \; t}}} \right)} + \ldots}} \\ {= {{\hat{\eta}}_{t + {k\; \Delta \; t}} + {{\hat{\phi}}_{1}e_{t}} + {{\hat{\phi}}_{2}e_{t - {\Delta \; t}}} + \ldots}} \end{matrix} & (15) \end{matrix}$

where “{circumflex over ( )}” is used for estimate and e_(t) is the Model 1 residual at t. Model 2, i.e., Eq. 14 provides estimation of the unknown parameters. While not wishing to be bound by the theory, we believe that the constraint given by Eq. 10b is very critical to achieving stable unbiased and accurate future predictions of BGC.

Applying this to the exemplary FBPC system in FIG. 1, y_(t) is provided in Eq. 14 and which theoretically should equal the set point (i.e., the desired BGC), Model 1 is provided in Eq. 8, and Model 2 is provided in Eq. 15.

In a further embodiment of the invention and as one of skill in the art would appreciate, BGC is affected by many different disturbances. The three most significant non-invasive disturbances fall into one of three categories: stress, consumed energy (i.e., meals), and activities. To adequately control BGC through the feedforward predictive control method, each of these disturbances must be taken into account. Announcements of each could be included in order to predict the impact of each as previously discussed. However, a benefit of the FBPC system and methods is that they overcome the difficulties previously encountered with the FFC system and methods which required modeling disturbances.

Method of Controlling Insulin

The process control of the invention can be employed in the monitoring and regulation of BGC for patients, particularly patients suffering from type 1 diabetes mellitus. The process control scheme of the invention can be applied to an insulin delivery system as shown in FIGS. 2A and 2B, which comprise a computing device 8 comprising a processor 10, a machine readable non-transitory media 12 which stores the FBPC control 14 of the invention, and a monitoring system 2.

Preferably, the method employs a monitoring system 2, which monitors the various inputs (disturbances). For example, the monitoring system can provide information regarding the activity inputs (X_(A) in FIGS. 3B and 3C), such as, skin temperature, near body temperature, and galvanic skin response. The monitoring system 2 can comprise one or more sensors capable of monitoring one of more of the inputs. In an aspect of the invention, the system can include a soft or virtual sensor monitoring, such as, heat flux. In an embodiment of the invention, the patient's basal and bolus are known from the insulin pump 2. For example, as bolus insulin can be administered by the insulin delivery system 4 and as the amount can be directed by the computing device 8, the amount of bolus insulin can be automatically logged. In a preferred embodiment, the methods are employed in an automatic insulin delivery system. An automatic insulin delivery system can be a system that monitors and instructs regarding insulin dosing, monitors and provides automatic insulin dosing, or an artificial pancreas.

As shown in FIGS. 2A and 2B, both the user input 16 and monitoring system 2 are in operable communication with a computing device 8 that contains a processor 10. The computing device 8 is also in operable communication with a machine readable non-transitory media 12. The machine readable non-transitory media 12 can store the FBPC control 14. The inputs 16, 2 are parametrized by the FBPC controller 14 to provide an output 6 that directs the insulin delivery system 4 to administer insulin. It is to be understood that wireless transceivers may be used to communicate information from sensors to the computing device 8. The sensors may be associated with wearable devices, mobile devices, or other computing devices.

FIGS. 2A and 2B show the computing device 8 is also in operable communication with an insulin delivery system 4. The insulin deliver system can comprise a wearable device, an artificial pancreas, or a combination thereof. The insulin delivery system 4 can comprise an apparatus suitable for administering insulin. Any suitable apparatus for administering insulin can be employed, including, but not limited to, an automatic insulin pump, a remotely controlled insulin pump, an IV, a catheter, or an artificial pancreas. Any suitable insulin pump capable of communication with the control system can be used. Preferably, the control method of the invention uses a remotely controlled insulin pump. Selection of an insulin pump can be based on preferences of the user for other features. The user inputs 16 and inputs provided by the monitoring system 2 are parametrized by the model 14 stored on the machine readable non-transitory media 12, which is in operable communication with the computing device 8. The FBPC control 14 provides an output to the computing device 8 which in turn provides an output 6 (as shown in FIG. 2A) to the insulin delivery system 4, which can direct the insulin delivery system 4 to adjust the insulin flow rate or quantity of insulin via the feedback predictive control scheme. In some embodiments of the invention, the output 6 can be provided in any communicable form and recorded, printed, displayed, etc. such that there is a record of the output 6. In an embodiment of the invention where the basal and bolus insulin are automatically monitored the output 6 can be relayed back into the machine readable non-transitory media 12.

Embodiments of the invention can reduce the variance in blood glucose concentration of a patient by the use of a more accurate future feedback error by predicting BGC more accurately due to a much short prediction horizon leading to better adjustments of insulin flow rate due under feedback predictive control. In embodiments of the invention, when compared with existing methods of controlling blood glucose concentration, the reduction in variance can be greater than 50%, preferably greater than 60%, more preferably greater than 70%, or most preferably greater than 80%.

EXAMPLES

Embodiments of the present invention are further defined in the following non-limiting Examples. It should be understood that these Examples, while indicating certain embodiments of the invention, are given by way of illustration only. From the above discussion and these Examples, one skilled in the art can ascertain the essential characteristics of this invention, and without departing from the spirit and scope thereof, can make various changes and modifications of the embodiments of the invention to adapt it to various usages and conditions. Thus, various modifications of the embodiments of the invention, in addition to those shown and described herein, will be apparent to those skilled in the art from the foregoing description. Such modifications are also intended to fall within the scope of the appended claims.

Example 1

Utilizing a continuous stirred tank reactor (CSTR) as a process model, the methods of the present invention were tested to verify Equations 1-15 and the FBPC algorithm exemplified in the block diagram of FIG. 1A. The FBPC algorithm was tested against an algorithm employing a model predictive control system as shown in FIG. 1B.

A multiple-input, single-output Wiener model as shown in FIG. 3A was used in modeling the CSTR process. Each input, X_(i), has its own linear dynamic block, G_(i), and each dynamic block has an intermediate, unobservable output V_(i), which represents the independent dynamic response of the corresponding input. The intermediate outputs are passed through a nonlinear static gain block, ƒ(V), to produce the final measured output, y.

Each dynamic block is a linear ordinary differential equation. In this Example, the manipulated variable, M, has first-order-plus-dead-time (FOPDT) dynamics shown in Eq. (16)

$\begin{matrix} {{{\tau_{1}\frac{{dv}_{1}(t)}{dt}} + {v_{1}(t)}} = {x_{1}\left( {t - \theta_{1}} \right)}} & (16) \end{matrix}$

where subscript ‘1’ denotes the manipulated variable, τ₁ is the time constant, and θ₁ is the dead time. Using a backward difference approximation, an approximate discrete-time form of Eq. (16) is given in Eq. (17)

v _(1,t)=δ_(1,1) v _(i,t−Δt)+(1−δ_(1,1))x _(1,t−k) ₁ _(Δt)  (17)

For k₁ steps into the future, Eq. (17) becomes Eq. (18):

$\begin{matrix} {{v_{1,{t + {k_{1}\Delta \; t}}} = {{\delta_{1,1}v_{1,{t + {{({k_{1} - 1})}\Delta \; t}}}} + {\left( {1 - \delta_{1,1}} \right)x_{1,t}}}}{where}} & (18) \\ {\delta_{1,1} = \frac{\tau_{1}}{\tau_{1} + {\Delta \; t}}} & (19) \end{matrix}$

The two other input variables, T_(ci) and C_(Ai), have second-order-plus-dead-time-plus-lead (SOPDTPL) form as shown in Eq. 20 below (which corresponds to Eq. 1 above):

$\begin{matrix} {{{\tau_{i}^{2}\frac{d^{2}{v_{i}(t)}}{{dt}^{\; 2}}} + {2\tau_{i}\zeta_{i}\frac{{dv}_{i}(t)}{dt}} + {v_{i}(t)}} = {{\tau_{ai}\frac{{dx}_{i}\left( {t - \theta} \right)}{dt}} + {x_{i}\left( {t - \theta_{i}} \right)}}} & (20) \end{matrix}$

where i=2 represents T_(ci), i=3 represents C_(Ai), τ_(i) is the time constant, ζ_(i) is the damping coefficient, τ_(ai) is the lead parameter, and θ_(i) is the dead time. Using a backward difference approximation applied to a sampling interval of Δt, an approximate discrete-time form of Eq. (20) is shown in Eq. (21) (corresponding to Eq. 2 above)

v _(i,t)=δ_(1,i) v _(i,t−Δt)+δ_(2,i) v _(i,t−2Δt)+ω_(1,i) x _(i,t−k) _(i) _(Δt)+ω_(2,i) x _(i,t−(k) _(i) _(+1)Δt)  (21)

where ω_(2,i)=1−δ_(1,i)−δ_(2,i)−ω_(1,i) to satisfy the unity gain constraint with

$\begin{matrix} {\delta_{1,i} = \frac{{2\tau_{i}^{2}} + {2\tau_{i}\zeta_{i}\Delta \; t}}{\tau_{i}^{2} + {2\tau_{i}\zeta_{i}\Delta \; t} + {\Delta \; t^{2}}}} & (22) \\ {\delta_{2,i} = \frac{- \tau_{i}^{2}}{\tau_{i}^{2} + {2\tau_{i}\zeta_{i}\Delta \; t} + {\Delta \; t^{2}}}} & (23) \\ {\omega_{1,i} = \frac{\left( {\tau_{ai} + {\Delta \; t}} \right)\Delta \; t}{\tau_{i}^{2} + {2\tau_{i}\zeta_{i}\Delta \; t} + {\Delta \; t^{2}}}} & (24) \end{matrix}$

Equations 22-24 correspond to Equations 3-5 above. After obtaining v_(i,t) for each input i, the modeled output value is determined by substituting these results into ƒ(V).

η_(t)=ƒ(V)=a ₀ +a ₁ v _(1,t) + . . . +a _(p) v _(p,t) +b ₁ v _(1,t) ² + . . . +b _(p) v _(p,t) ² +c _(1,2) v _(1,t) v _(2,t) + . . . +c _(p−1,p) v _(p−1,t) v _(p,t)  (25)

Modification of Eq. (25) for predicting k₁ time steps into the future with p=3 gives

η_(t+k) ₁ _(Δt) =a ₀ +a ₁ v _(1,t+k) ₁ _(Δt) + . . . +a ₃ v _(3,t+k) ₁ _(Δt)+

b ₁ v _(1,t+k) ₁ _(Δt) ² + . . . +b ₃ v _(3,t+k) ₁ _(Δt) ² +c _(1,2) v _(1,t+k) ₁ _(Δt) v _(2,t+k) ₁ _(Δt) + . . . +c _(2,3) v _(2,t+k) ₁ _(Δt) v _(3,t+k) ₁ _(Δt)  (26)

where a_(i), b_(i), and c_(i,j), denote the linear, quadratic and interaction parameters for i=1, 2, 3 and j=2 and 3, and

v _(2,t+k) ₁ _(Δt)=δ_(1,2) v _(2,t+(k) ₁ _(−1)Δt)+δ_(2,2) v _(2,t+(k) ₁ _(−2)Δt)+ω_(1,2) x _(2,t+(k) ₁ _(−k) ₂ _()Δt)+ω_(2,2) x _(2,t+(k) ₁ _(−k) ₂ _(−1)Δt)  (27)

v _(3,t+k) ₁ _(Δt)=δ_(1,3) v _(3,t+(k) ₁ _(−1)Δt)+δ_(2,3) v _(3,t+(k) ₁ _(−2)Δt)+ω_(1,3) x _(3,t+(k) ₁ _(−k) ₃ _()Δt)+ω_(2,3) x _(3,t+(k) ₁ _(−k) ₃ _(−1)Δt)  (28)

In this Example, all measured inputs are assumed to be known over the prediction horizon.

To model unmeasured disturbances and bias, a noise model structure was applied to Eq. (25) as follows:

$\begin{matrix} {\mspace{79mu} {{y_{t} = {\eta_{t} + N_{t}}},\mspace{79mu} {where}}} & (26) \\ {N_{t} = {{\frac{\theta_{q}(B)}{\varphi_{p}(B)}a_{t}} = {{\frac{1 - {\theta_{1}B} - {\theta_{2}B^{2}} - \ldots - {\theta_{q}B^{q}}}{1 - {\varphi_{1}B} - {\varphi_{2}B^{2}} - \ldots - {\varphi_{p}B^{p}}}a_{t}} = {\frac{a_{t}}{\Phi (B)} = {\frac{a_{t}}{\left\lbrack {1 - {\phi_{1}B} - {\phi_{2}B^{2}} - \ldots}\mspace{11mu} \right\rbrack}{\forall t}}}}}} & (27) \\ {\mspace{79mu} {a_{t}\overset{indep}{\sim}{N\left( {0,\sigma^{2}} \right)}}} & (28) \end{matrix}$

and y_(t) is the measured tank temperature, T_(m), at t and B^(γ)x_(t)=X_(t−γΔt).

Then

$\begin{matrix} {{y_{t}\left\lfloor {1 - {\phi_{1}B} - {\phi_{2}B^{2}} - \ldots}\mspace{11mu} \right\rfloor} = \left. {{\eta_{t}\left\lfloor {1 - {\phi_{1}B} - {\phi_{2}B^{2}} - \ldots}\mspace{11mu} \right\rfloor} + a_{t}}\Rightarrow \right.} & (29) \\ {y_{t} = {\eta_{t} + {\phi_{1}\left( {y_{t - {\Delta \; t}} - \eta_{t - {\Delta \; t}}} \right)} + {\phi_{2}\left( {y_{t - {2\Delta \; t}} - \eta_{t - {2\Delta \; t}}} \right)} + \ldots + a_{t}}} & (30) \end{matrix}$

Modification of Eq. (30) for predicting k₁ time steps into the future gives

$\begin{matrix} \begin{matrix} {{\hat{y}}_{t + {{k\;}_{1}\Delta \; t}} = {{\hat{\eta}}_{t + {k_{1}\; \Delta \; t}} + {{\hat{\phi}}_{1}\left( {y_{t} - {\hat{\eta}}_{t}} \right)} + {{\hat{\phi}}_{2}\left( {{\hat{y}}_{t - {\Delta \; t}} - {\hat{\eta}}_{t - {\Delta \; t}}} \right)} + \ldots}} \\ {= {{\hat{\eta}}_{t + {{k\;}_{1}\Delta \; t}} + {{\hat{\phi}}_{1}e_{t}} + {{\hat{\phi}}_{2}e_{t - {\Delta \; t}}} + \ldots}} \end{matrix} & (31) \end{matrix}$

where “{circumflex over ( )}” is used for estimate. {circumflex over (η)}_(t+k) ₁ _(Δt) is called Model 1 and ŷ_(t+k) ₁ _(Δt) is called Model 2. In this Example, since the input values, as well as the current and past values of the controlled variable, are known, Eq. 31 has all the necessary information to predict the value of the controlled variable at t+k₁Δt (i.e., t+θ₁) time in the future.

The block diagrams of the two virtual sensor predictive control algorithms are shown in FIGS. 1 & 2. A description of the transfer functions in the blocks is given in Table 1.

TABLE 1 Transfer Function Description G_(c) Controller G_(v) Valve G_(p) Process G_(Li) Load i G_(m) Transmitter for the controlled variable K_(m) Transmitter gain for the controlled variable

For each of the systems, a change in the manipulated variable at t affects the controlled variable k₁Δt time into the future. As shown in FIG. 1, for FBPC the controller is a common PID controller, and thus can have a three-dimensional tuning space. The FBE requires the predicted value of the controlled variable k₁Δt time into the future, the time when current change in the manipulated variable affects the output. In contrast to FBPC, MPC has one tuning parameter, J, representing the number of time steps beyond t+k₁Δt when the multiple-input model of the controlled variable (i.e., ŷ_(t+(k) ₁ _(+J)Δt)) equals the set point for the next change. Thus, MPC requires prediction farther into the future than FBPC, and hence, its FBE, will likely be less accurate. Also, as shown in FIG. 2, the MPC control law uses the inverse of a one-input controlled variable model (Ĝ_(VP)), to determine the next control move. Hence, it is critical that this model has an accurate causative relationship, which can be difficult to achieve in practice.

The CSTR process used in this Example is illustrated in FIG. 4 and described by Eqs. 32-36. The definitions and values of the variables are outlined in Table 2.

$\begin{matrix} {\frac{{dC}_{A}}{dt} = {{\frac{F}{V}\left( {C_{Ai} - C_{A}} \right)} - {kC}_{A}^{2}}} & (32) \\ {\frac{{dT}_{T}}{dt} = {{\frac{F}{V}\left( {T_{T} - T_{T}} \right)} - {\frac{\Delta \; H_{R}}{\rho \; C_{p}}{kC}_{A}^{2}} - {\frac{UA}{\rho \; {VC}_{p}}\left( {T_{T} - T_{C}} \right)}}} & (33) \\ {\frac{{dT}_{C}}{dt} = {{\frac{UA}{V_{C}\rho_{C}C_{pC}}\left( {T_{T} - T_{C}} \right)} - {\frac{F}{V_{C}}\left( {T_{C} - T_{Ci}} \right)}}} & (34) \\ {k = {k_{0}\mspace{11mu} {\exp \left( \frac{- E}{R\left( {T_{T} + 273.16} \right)} \right)}}} & (35) \\ {F_{C} = {F_{C\mspace{11mu} \max} \times \alpha^{({- M})}}} & (36) \end{matrix}$

TABLE 2 Variable Definition SS value (units) A Heat transfer area 5.40 (m²) α Control valve 50 (none) rangeability parameter C_(A) Concentration of 1.0302 (kgmol/m³) species A in reactor C_(Ai) Concentration of species 2.88 (kgmol/m³) A in inlet stream c_(p) Heat capacity of feed and 1.815 × 10⁵ (J/kgmol-° C.) product streams c_(pc) Heat capacity of coolant 4184 (J/kg-° C.) ΔH_(R) Heat of reaction −9.86 × 10⁷ (J/kgmol) E Activation energy 1.182 × 10⁷ (J/kgmol) F Feed flow rate 0.45 (m³/s) F_(C) Coolant flow rate 0.44 (m³/s) F_(Cmax) Maximum flow rate of 1.2 (m³/s) coolant through control valve K Reaction rate constant 0.09 (m³/s-kgmol) k_(o) Arrhenius frequency 0.0744 (m³/s-kgmol) parameter M Input signal to the valve 0.26 (none) subject to 0 ≤ M ≤ 1 R Gas law constant 8314.39 (J/kgmol-K) ρ Density of reactor 19.2 (kgmol/m³) contents ρ_(c) Density of coolant 1000 (kg/m³) T_(c) Coolant temperature 50.48 (° C.) in the jacket T_(Ci) Coolant inlet 27 (° C.) temperature T Reactor temperature 88 (° C.) T_(m) Measured reactor 88 (° C.) temperature U Overall heat transfer 2.13 × 10⁵ (J/s-m²-° C.) coefficient V_(C) Cooling jacket volume 1.82 (m³) V CSTR volume 7.08 (m³)

In this Example, three input variables are used in the model. They are: 1) the input signal to the valve, M; 2) the temperature of the entering coolant, T_(ci); and 3) the concentration of component A in the inlet stream, C_(Ai). There is one unmeasured disturbance, the feed flow rate, F. There is one output for this Example, the measured tank temperature, T_(m). The CSTR is altered to include dead time, θ_(t), in T_(ci), C_(Ai), and the input signal to valve, M, specified as θ_(i)=k_(i)Δt, where k_(i) is a positive integer, i=1, 2, or 3, for M, T_(ci), and C_(Ai), respectively. The dead time forM is the largest.

The training input sequences for M, T_(ci) and C_(Ai) are shown in FIG. 5. The fits for this training data are shown in FIGS. 6 and 7. FIG. 6 shows the excellent fit of Model 2, the model that fits the inputs and the serial correlation structure of the noise, the unmeasured disturbances, and the model bias, as given by Eq. (27). This is the model used in the control algorithms. The fit shown in FIG. 7 is a fit of the input-only Model 1. Thus, it does not fit the unmeasured disturbance, and as a result, this fit is very poor, given the large changes in the feed rate, F.

In this Example, k₁Δt, k₂Δt and k₃Δt are 10 seconds, 5 seconds, and 0 seconds with Δt=0.1 seconds, thus the predictions for each control method are at least 100 time steps into the future. FIG. 8A shows the response of the measured reactor temp, T_(m), the variable to be controlled, when the input signal to the valve, M, is fixed at its initial value (as shown in FIG. 8B). The other input variables are changed as shown in FIG. 8B as if for the control run.

The objective is to minimize the variation of T_(m) around the set point by manipulating the input signal to the valve, M, subject to 0≤M≤1. To correct for bias, the common practice for MPC applications is to use a bias correction that is the latest measurement of T_(m) minus its corresponding predicted value. When the MPC estimator was used in the bias correction, the results were poor. This is shown as “MPC-BC” in Table 3. In addition, when these differences as the residuals were used for Model 2 in the MPC algorithm, control was also poorer than FBPC. The best results for MPC were obtained using Model 2 to correct for bias. These control results along with the ones for FBPC are given in Table 3.

TABLE 3 Control Stdev % Increase Method J (° C.) from FBPC FBPC NA 0.41 0 MPC 10 0.49 18 15 0.49 18 20 0.49 18 30 0.48 17 MPC—BC 30 0.96 132

As shown in Table 3, FBPC has the smallest Stdev. MPC for J=10, 15, 20, and 30 are all similar and about 20% worse than FBPC. Note that for these values of J, MPC is only predicting 1 to 3 seconds farther into the future. Thus, for such a short distance into the future and the high accuracy of the fitted model, the improvement in FBPC is not expected to be too much greater than MPC, as the results indicate. T_(m) and M results for FBPC are shown in FIGS. 9A and 9B, respectively. The T_(m) and M results for MPC where J=10 and 15, are shown in FIGS. 9C-9F. In going from J=10 to 15, M for MPC is more aggressive than FBPC. Thus, even with more aggressiveness in the manipulated variable, MPC control performance is worse than FBPC. The MPC case with bias correction, representing the current practice, denoted by ‘MPC-BC’ in Table 3, is 132% worse than FBPC. Thus, Model 2 helped the control of MPC substantially over the bias correction method but MPC results are still worse than FBPC.

Example 2

In this Example, the present invention's modeling approach is shown to accurately predict BGC 30 and 60 minutes into the future for Type 1 diabetes subjects, which is necessary to obtain an accurate prediction of FBE to effectively manipulate insulin flow rate at the present time for optimal control of BGC. This Example is an open-loop modeling study, and not a closed-loop control study.

For this Example, the process model utilizes a Wiener model with a semi-coupled network for BIC and BGC as illustrated in FIG. 3B. The semi-coupled network allows modeling of blood insulin infusion rate and glucose interaction. BIC is not measured, and therefore the approach uses a “pseudo” BIC variable that allows the use of a modeling structure that is hypothesized to provide characteristics for BGC that are more physiologically correct than a Wiener-only modeling technique.

This Example employs the two-step modeling approach outlined in Example 1: applying a noise model structure (Model 2) to model unmeasured disturbances and bias in the predictive model (Model 1). Model 1 uses the network illustrated in FIG. 3B. As shown all inputs, X_(i), pass through a linear dynamic block to produce the unobservable dynamic output variables, V_(i). Note that i=A_(i) . . . A_(p) for p activity inputs, C for carbohydrates, F for fats, P for proteins, and I for insulin. The dynamic blocks for the inputs have SOPDTPL for as shown in Example 1 Eqs. (20)-(24). A dynamic mass balance on the G_(FI) block in FIG. 3B, which represents the BGC due to food and insulin only, is shown in Eq. (37)

$\begin{matrix} {\frac{{dV}_{FI}(t)}{dt} = {{a_{C}{v_{C}(t)}} + {a_{F}{v_{F}(t)}} + {a_{P}{v_{P}(t)}} - {a_{FI}{I(t)}{V_{FI}(t)}}}} & (37) \end{matrix}$

where a_(C), a_(F), and a_(FI), are estimated model parameters, I(t) is the unmeasured BIC at time t, and V_(FI)(t) is the BGC due to food and insulin changes only. Similarly, a dynamic mass balance on the G_(BI) block, which represents BIC is shown in Eq. (38)

$\begin{matrix} {\frac{{dI}(t)}{dt} = {{a_{I}{v_{I}(t)}} - {a_{BI}{I(t)}}}} & (38) \end{matrix}$

where a₁ and a_(BI) are estimated model parameters. The function ƒ(V) is a first-order linear regression structure given in discrete form in Eq. (39)

ƒ(V)=η_(t)=ε_(i) =a ₀ +V _(FI,t) +a _(A1) v _(A) _(t) _(,t) + . . . +a _(Ap) v _(A) _(p) _(,t)  (39)

where a₀ and c_(Ai)'s are static estimated parameters.

Unmeasured disturbances and bias, Model 2, is as described mathematically in Example 1, Eqs. (26)-(31), subject to one constraint as described in Eqs. (40) and (41) (corresponding to equations 10a and 10b above), which is critical to achieving stable unbiased and accurate future prediction of BGC.

$\begin{matrix} {N_{t} = {{\frac{\theta_{q}(B)}{\varphi_{p}(B)}a_{t}} = {{\frac{1 - {\theta_{1}B} - {\theta_{2}B^{2}} - \ldots - {\theta_{q}B^{q}}}{1 - {\varphi_{1}B} - {\varphi_{2}B^{2}} - \ldots - {\varphi_{p}B^{p}}}a_{t}} = {\frac{a_{i}}{\Phi (B)} = {\frac{a_{t}}{\left\lbrack {1 - {\phi_{1}B} - {\phi_{2}B^{2}} - \ldots}\mspace{11mu} \right\rbrack}{\forall t}}}}}} & (40) \\ {\mspace{79mu} {{{subject}\mspace{14mu} {to}\mspace{14mu} {\sum\phi_{i}}} = 1}} & (41) \end{matrix}$

In this Example, the model was fit to eleven cases of BGC subject data split into sets of training, validation, and testing. The model parameters are estimated using the training data, and the validation set is used to guard against over-fitting, as the final set of parameters must maximize the fit in this set. Since these two sets influence the parameter estimation process, the testing set is used as a final check on model fit, as the data in this set has no influence on the values of the estimated parameters. This Example applies Model 2 under the following conditions for the training data sets:

$\begin{matrix} {{\min\limits_{{\hat{\phi}}_{i}}\mspace{11mu} {SSE}} = {\sum\limits_{{over}{training}{time}}e_{t}^{2}}} & (42) \\ {{{subject}\mspace{14mu} {to}\mspace{14mu} {\sum\limits_{i = 1}^{4}{\hat{\phi}}_{i}}} = 1} & (43) \end{matrix}$

Thus, four {circumflex over (φ)} were used in each case and convergence criterion is given in Eq. (42) and subject to Eq. (43). In this Example using Model 2, all the data not used to determine the {circumflex over (φ)}_(i)'s is testing data.

The results for this Example are reported using three metrics. The first is the correlation coefficient between the fitted and measured BGC, r_(fit) and is given by Eq. (44):

$\begin{matrix} {r_{fit} = \frac{{\sum\limits_{i = 1}^{n_{t}}{y_{i}{\hat{y}}_{i}}} - \frac{\left( {\sum\limits_{i = 1}^{n_{t}}y_{i}} \right)\left( {\sum\limits_{i = 1}^{n_{t}}{\hat{y}}_{i}} \right)}{n_{t}}}{\sqrt{{\sum\limits_{i = 1}^{n_{t}}y_{i}^{2}} - \frac{\left( {\sum\limits_{i = 1}^{n_{t}}y_{i}} \right)^{2}}{n_{t}}}\sqrt{{\sum\limits_{i = 1}^{n_{t}}{\hat{y}}_{i}^{2}} - \frac{\left( {\sum\limits_{i = 1}^{n_{t}}{\hat{y}}_{i}} \right)^{2}}{n_{t}}}}} & (44) \end{matrix}$

where n_(t) is the number of pairs of values in the set, and y_(i) and ŷ_(i) are the i^(th) measured and fitted values in the set, respectively. The closer this statistic is to 1.0, the better the model fits the measured response data without consideration of model bias.

The second statistic is the averaged difference, AD, which gives a measure of model bias. As shown in Eq. (45), it is the average of deviation for the differences of y_(i) and y_(i).

$\begin{matrix} {{AD} = {\frac{1}{n_{t}}{\sum\limits_{i = 1}^{n_{t}}\left( {y_{i} - {\hat{y}}_{i}} \right)}}} & (45) \end{matrix}$

The final statistic is the averaged absolute error, AAD, and is a measure of the average closeness of y_(i) and ŷ_(i), as shown in Eq. (46).

$\begin{matrix} {{AAD} = {\frac{1}{n_{t}}{\sum\limits_{i = 1}^{n_{t}}{{y_{i} - {\hat{y}}_{i}}}}}} & (46) \end{matrix}$

The 13 variable input set is given in Table 4 and consists of three food nutrients, seven activity variables and two for insulin infusion. The activity variables were collected using the SenseWear® Pro3 Body Monitoring System shown in FIG. 14, which is worn on the triceps of the subject's arm. The SenseWear® armband utilizes pattern detection algorithms that employ physiologic signals from a unique combination of sensors to generate values for 20 activity variables. The armband collects data using a two-axis accelerometer and four sensors that are used to determine heat flux, skin temperature, near body temperature, and galvanic skin response. The two-axis accelerometer provides information about body position and tracks upper arm movement. The heat flux sensor calculates the amount of heat being dissipated from the body by measuring the amount of heat lost along a thermally-conductive path between the skin and a vent on the side of the armband. Skin temperature and near-body temperature are measured by sensitive thermistors, and galvanic skin response is measured via the conductivity of the subject's skin as it varies due to physical and emotional stimuli. Measurements were taken at five-minute intervals. The data collection was conducted over the course of a year by an experienced graduate student.

TABLE 4 Food Inputs Activity Inputs Circadian Rhythm Insulin Inputs 1 = Carbohydrates 4 = Transverse 11 = Time of Day 12 = Bolus acceleration (peaks) (TOD) 2 = Fats 5 = Heat Flux (avg) 13 = Basal 3 = Proteins 6 = Longitudinal acceleration (avg) 7 = Near body temperature 8 = Transverse acceleration (MAD) 9 = GSR (avg) 10 = Energy expenditure

The results of the predictive modeling approach 30 and 60 minutes into the future are outlined in Tables 5 and 6, respectively, with the data split into three sets comprising of one week of training, four days of validation, and three days of testing. Tables 7 and 8 show results for modeling 30 and 60 minutes into the future, respectively, with the data split into two sets: one week of training and one week of validation.

TABLE 5 7 Day Training 4 Day Validation 3 Day Test Subject AD AAD rfit AD AAD rfit AD AAD rfit 1 −0.4 16.2 0.95 −0.3 19.5 0.93 0.7 15.8 0.96 2 0.3 20.6 0.96 −1.6 22.0 0.94 −0.5 26.7 0.94 3 0.8 24.0 0.92 −0.6 25.4 0.89 0.5 22.6 0.90 4 −0.1 17.4 0.87 −0.3 23.4 0.90 0.1 18.0 0.89 5 0.1 26.2 0.92 1.2 29.6 0.94 −2.2 26.4 0.92 6 0.0 17.2 0.95 0.6 14.9 0.97 −0.1 14.8 0.93 7 −0.1 21.3 0.93 −0.9 19.6 0.87 0.9 19.6 0.95 8 −0.1 15.3 0.90 0.1 22.7 0.83 1.8 18.3 0.90 9 0.4 18.7 0.95 0.2 23.5 0.91 0.0 18.1 0.96 10 −0.3 18.4 0.95 −0.6 17.0 0.95 0.6 25.5 0.92 11 −0.3 17.9 0.93 1.4 20.3 0.93 −0.7 19.4 0.95 Average 0.0 19.4 0.93 −0.1 21.6 0.91 0.1 20.5 0.93

TABLE 6 7 Day Training 4 Day Validation 3 Day Test Subject AD AAD rfit AD AAD rfit AD AAD rfit 1 −1.2 29.0 0.84 −1.6 31.3 0.85 1.4 26.4 0.91 2 0.9 33.5 0.91 −5.1 34.8 0.86 −1.9 42.2 0.86 3 1.3 37.2 0.82 −1.1 42.6 0.72 0.9 37.9 0.77 4 −0.3 25.5 0.73 −0.4 38.8 0.70 −0.1 29.4 0.76 5 0.5 41.6 0.82 2.0 48.4 0.85 −3.6 36.3 0.87 6 0.7 31.1 0.86 0.8 25.0 0.91 −0.3 25.1 0.80 7 −0.5 37.0 0.83 −2.4 37.2 0.70 2.1 37.2 0.85 8 −0.1 23.9 0.77 0.1 29.9 0.74 3.3 28.0 0.80 9 0.6 30.6 0.87 1.0 37.2 0.82 0.7 31.2 0.88 10 −0.7 30.1 0.86 −0.6 30.1 0.85 1.2 38.1 0.85 11 −1.0 26.1 0.86 2.5 29.3 0.85 −1.3 29.6 0.90 Average 0.0 31.4 0.83 −0.4 35.0 0.80 0.2 32.9 0.84

TABLE 7 Training Validation Subject AD AAD rfit AD AAD rfit 1 −0.4 15.5 0.95 0.5 18.1 0.94 2 0.5 25.0 0.95 1.7 24.3 0.94 3 0.5 23.2 0.93 0.5 22.2 0.91 4 0.0 17.5 0.88 −0.3 20.4 0.90 5 −0.3 25.3 0.93 0.1 25.9 0.95 6 0.3 17.5 0.95 0.3 14.6 0.96 7 −0.1 20.0 0.94 0.4 18.8 0.93 8 −0.1 15.3 0.90 0.8 20.6 0.87 9 −0.1 20.9 0.94 0.4 24.3 0.92 10 −0.3 18.6 0.94 0.1 22.5 0.93 11 −0.3 17.4 0.93 0.5 17.9 0.95 Average 0.0 19.7 0.93 0.5 20.9 0.93

TABLE 8 Training Validation Subject AD AAD rfit AD AAD rfit 1 −1.2 28.5 0.86 0.4 32.8 0.85 2 1.1 37.8 0.90 2.7 36.4 0.87 3 0.7 36.3 0.84 1.2 35.2 0.79 4 −0.1 25.2 0.73 −0.9 33.0 0.74 5 −0.2 42.1 0.82 0.4 43.8 0.85 6 1.3 31.9 0.85 0.5 24.6 0.88 7 −0.3 33.8 0.84 1.0 30.6 0.85 8 −0.1 23.6 0.76 1.5 28.5 0.78 9 −0.5 35.2 0.85 1.5 41.5 0.8 10 −0.5 30.3 0.85 0.7 35.9 0.84 11 −0.8 25.3 0.86 1.0 27.4 0.89 Average −0.1 31.8 0.83 0.9 33.6 0.83

The results of the predictive modeling method as outlined in Tables 5-8 indicate very good to excellent outcomes. The models have little to negligible bias as shown by the very low AD values. For the 30-minute predictive models, AAD averages just over 20 mg/dL. For the 60-minute models, the AAD increases about 12 mg/dL. For the 30-minute prediction study, the average r_(fit) is about 0.93, which is excellent, with highs of 0.96. For the 60-minute study, the drop is about 0.1 on the averages to about 0.83. However, there are some with r_(fit) of 0.9 or better.

Example 3

In this Example, a diabetes simulator developed by UVA/Padova, and approved by the FDA, was used to test the present invention. This simulator provides a non-human methodology to test different BGC control algorithms and is considered a valid substitute for animal trials in the early stages for development of an artificial pancreas. This diabetes simulator was used to generate subject data for a comparison of MPC and FBPC with all things being equal except for the differences in the control algorithms. The fitted predictive model used by both methods is the same and predicts G accurately, thus any difference in control performance is due to differences in the control algorithms.

The FBPC scheme is illustrated in FIG. 13 and is described in FIG. 2. It should be noted that in Example 1 the FBPC method was used in a CSTR system and in this Example the inputs and outputs are specific to controlling BGC by varying insulin flow rate. The present system utilizes a Wiener Model as shown in FIG. 1 and described mathematically in Example 1. In the general context of diabetes, the inputs, X for i=1 . . . p, are the measured noninvasive variables such as meal components, physical activity, and emotional stress. The output, y, is BGC. In the diabetes simulator used in this Example, there are no unmeasured disturbances and only two inputs: the amount of carbohydrates in each meal, and IFR.

This Example uses 30 virtual patients, the characteristics of which are summarized in Table 9. To minimize the effect of measurement delay, both the insulin pump and BGC sensor were selected as intravenous (IV) type. As for the scenario tested in this Example, three meals occurred each day at 7 am, 12 pm and 8 pm, and all meals were 50 grams of carbohydrate. To examine the robustness of the proposed FBPC method, a longer than common period of a five-day closed loop control run was utilized. At the start of each trial, a two-day “run-in” period was used to bring the simulator to a consistent day-to-day behavior. Therefore, in total, a seven-day run was used for each subject.

TABLE 9 Body Weight Fasting BGC Subject Age (kg) (mg/dL) 1 16 36.16 115 2 18 50.88 123 3 19 54.49 127 4 17 41.47 123 5 19 59.34 118 6 17 52.58 119 7 19 48.58 133 8 16 51.79 112 9 16 67.01 120 10 19 49.36 117 11 3 27.79 123 12 12 32.61 124 13 5 35.20 111 14 5 28.32 124 15 3 36.99 121 16 2 31.62 118 17 2 25.48 122 18 8 28.86 123 19 4 23.01 125 20 7 36.15 133 21 32 79.80 122 22 22 80.34 123 23 42 70.87 119 24 24 67.48 119 25 47 67.11 120 26 23 73.01 116 27 47 46.06 115 28 56 98.66 130 29 24 67.78 114 30 31 80.79 113

Both control algorithms were supplemented by basal insulin infusion, at the rate of U/hr, that is customized to each of the test subjects as given in Table 9. The set point for each subject was set to their fasting BGC before the start of the trial.

The manipulated input variable, insulin feed rate (IFR), is atypical because it can only decrease BGC. Thus, when a control system calculated the value for bolus insulin to be negative, bolus IFR was necessarily set to zero. In addition, for FBPC, to avoid hypoglycemia, the implemented protocol was to turn bolus IFR off one minute after any meal and then turn it back on one hour before the meal.

A common benchmark used in this area is the percent of values in a range, in this Example, between 70-180 mg/dL as suggested by American Diabetes Association. Values below 70 mg/dL are commonly considered as the hypoglycemic region and above 180 mg/dL is commonly considered as the hyperglycemic region. Thus, the effectiveness a BGC control algorithm is commonly reported by the percentage of time spent in 70-180 mg/dL range. Additionally, other metrics can describe different aspects of the effectiveness such as the spread of BGC about its mean, or its standard deviation (Stdev). In this Example, both metrics are reported.

Three measures of performance are used in this Example for evaluating the fit of the model. The first is the correlation coefficients between the fitted and measured BGC (r_(fit)) and given by Eq. (44):

$\begin{matrix} {r_{fit} = \frac{{\sum\limits_{i = 1}^{n_{t}}{y_{i}{\hat{y}}_{i}}} - \frac{\left( {\sum\limits_{i = 1}^{n_{t}}y_{i}} \right)\left( {\sum\limits_{i = 1}^{n_{t}}{\hat{y}}_{i}} \right)}{n_{t}}}{\sqrt{{\sum\limits_{i = 1}^{n_{t}}y_{i}^{2}} - \frac{\left( {\sum\limits_{i = 1}^{n_{t}}y_{i}} \right)^{2}}{n_{t}}}\sqrt{{\sum\limits_{i = 1}^{n_{t}}{\hat{y}}_{i}^{2}} - \frac{\left( {\sum\limits_{i = 1}^{n_{t}}{\hat{y}}_{i}} \right)^{2}}{n_{t}}}}} & (44) \end{matrix}$

where n_(t) is the number of pairs of values in the set, and y_(i) and ŷ_(i) are the i^(th) measured and fitted values in the set, respectively. The closer this statistic is to 1.0, the better the model fits the measured data without consideration of model bias. The second statistic is the averaged error (AE), which gives a measure of model bias. As shown in the equation below, it is the average of y_(i)−ŷ_(i).

$\begin{matrix} {{AE} = {\frac{1}{n_{t}}{\sum\limits_{i = 1}^{n_{t}}\left( {y_{i} - {\hat{y}}_{i}} \right)}}} & (45) \end{matrix}$

The final statistic is the averaged absolute error (AAE), which is a measure of the average closeness of y_(i) to ŷ_(i), as shown in Eq. (46).

$\begin{matrix} {{AAE} = {\frac{1}{n_{t}}{\sum\limits_{i = 1}^{n_{t}}{{y_{i} - {\hat{y}}_{i}}}}}} & (46) \end{matrix}$

Model development followed the cross-validation protocol of Example 1 as described in Eq. (30). The size and type of data sets for model identification are the same as in Eq. (26). The controllers were tuned to give the tightest control of glucose levels. Table 10 gives a summary of the modeling results in this Example for each subject. As shown, r_(fit) is at, or nearly at, its maximum value of one (1) for all 30 subjects. Similarly, the AE and AAE results are also excellent averaging 0.22 and 1.7 mg/dL, respectively. To demonstrate the excellent accuracy of the model n time steps into the future under automatic control, results are shown in FIG. 10 for Subject 4 with r_(fit)=0.99, AE=0.068, and AAE=1.897.

FIGS. 11A-E show graphical control run results for both algorithms (MPC and FBPC) for the first five (5) subjects and FIGS. 12A-E show insulin feed rates (IFR) for these subjects during the control runs. The control run results for all the subjects are outlined in Tables 11A and 11B. Table 12 gives the tuning parameters for the first five (5) subjects.

TABLE 10 AE AAE r_(fit) (mg/dL) (mg/dL) τ₁ τ₂ Sub- Train- Valida- Test- Test- Test- Carb Insulin ject ing tion ing ing ing (min) (min) 1 0.999 0.999 0.999 0.275 1.694 574.9 804.4 2 0.997 0.996 0.997 0.433 2.161 654.9 916.5 3 0.999 0.999 0.999 0.219 1.218 669.0 871.5 4 0.999 0.999 0.999 0.068 1.897 796.8 1131.9 5 1.000 0.999 0.999 0.088 1.454 1015.6 1078.1 6 0.999 0.999 1.000 0.271 0.916 715.3 1145.6 7 1.000 1.000 1.000 0.123 1.110 571.0 1078.4 8 0.999 0.999 0.999 0.510 1.502 813.8 800.0 9 0.999 0.999 0.999 0.226 1.456 851.0 1048.9 10 0.998 0.997 0.998 0.242 1.692 864.5 1112.7 11 0.999 0.999 0.999 0.878 2.890 736.1 818.0 12 0.999 0.999 0.999 0.589 2.243 613.1 879.6 13 0.999 0.999 0.999 0.574 2.361 723.3 943.4 14 0.999 0.999 0.999 1.790 4.066 905.6 994.6 15 0.996 0.995 0.997 0.437 3.335 540.7 772.5 16 0.998 0.998 0.999 0.616 3.359 680.2 814.9 17 0.998 0.998 0.999 1.446 5.032 782.8 806.2 18 0.998 0.998 0.999 0.426 3.532 822.4 1032.4 19 0.999 0.999 0.999 1.416 3.982 927.8 914.8 20 0.999 0.998 0.999 0.625 2.422 779.0 916.3 21 0.997 0.997 0.998 0.129 1.271 630.7 821.2 22 0.996 0.996 0.997 0.043 0.871 447.9 578.2 23 0.999 0.999 0.999 0.231 1.061 779.0 883.6 24 0.999 0.998 0.999 0.093 1.076 681.8 1053.9 25 0.999 0.999 0.999 0.148 1.190 798.9 1063.5 26 0.998 0.997 0.998 0.049 1.151 480.0 800.0 27 0.999 0.999 0.999 0.705 2.319 867.2 1045.1 28 0.998 0.998 0.998 0.264 1.521 664.9 902.8 29 0.999 0.998 0.999 0.115 1.138 549.5 813.2 30 0.998 0.997 0.998 0.111 1.028 580.8 780.0 Mean 0.999 0.998 0.999 0.438 2.032 717.3 920.7

TABLE 11A FBPC MPC 70-180 70-180 Set Mean Stdev mg/dL Mean Stdev mg/dL Subject Point (mg/dL) (mg/dL) (%) (mg/dL) (mg/dL) (%) 1 115 117.9 34.2 94.7 108.2 42.4 69.8 2 123 123.4 22.5 97.8 110.8 32.7 94.2 3 127 127.3 20.6 100.0 114.8 23.6 100.0 4 123 128.1 33.4 96.0 130.0 41.3 84.2 5 118 118.1 13.1 100.0 108.1 21.9 100.0 6 119 119.4 20.8 100.0 106.9 25.5 100.0 7 133 134.4 42.3 78.0 104.3 48.0 60.7 8 112 112.0 20.6 99.9 103.1 24.9 93.8 9 120 119.9 17.4 100.0 115.5 23.0 100.0 10 117 116.9 16.6 100.0 109.3 25.6 98.7 11 123 123.1 35.0 92.4 121.9 39.0 87.6 12 124 125.9 40.6 83.7 112.2 47.9 65.6 13 111 127.4 29.4 96.1 107.2 33.5 85.0 14 124 123.7 22.2 97.7 105.9 28.7 89.0 15 121 121.6 32.9 93.6 117.9 41.9 78.7 16 118 120.7 21.0 100.0 123.2 34.8 88.5 17 122 122.6 26.5 93.0 126.9 29.6 91.8 18 123 122.0 36.1 92.1 97.9 40.8 65.0 19 125 124.8 13.7 100.0 120.8 17.9 100.0 20 133 133.1 20.4 98.4 123.8 26.7 97.0 21 122 122.0 13.4 100.0 121.6 19.1 100.0 22 123 123.6 9.6 100.0 123.5 15.9 100.0 23 119 121.5 11.1 99.6 119.9 15.8 99.7 24 119 119.1 14.1 100.0 115.7 18.6 100.0 25 120 120.0 11.9 100.0 120.6 16.4 100.0 26 116 121.1 16.9 100.0 121.4 24.2 99.1 27 115 115.0 14.3 100.0 112.1 22.5 100.0 28 130 130.0 14.5 100.0 128.3 19.5 100.0 29 114 112.9 15.5 100.0 113.5 22.1 100.0 30 113 112.7 13.5 100.0 109.0 16.2 100.0 Average 122.0 21.8 97.1 115.1 28.0 91.6

TABLE 11B MPC MPC Stdev > Stdev − FPC FBPC Set Stdev Stdev Subject Point (%) (mg/dL) 1 115 24.0 8.2 2 123 45.3 10.2 3 127 14.6 3.0 4 123 23.7 7.9 5 118 67.2 8.8 6 119 22.6 4.7 7 133 13.5 5.7 8 112 20.9 4.3 9 120 32.2 5.6 10 117 54.2 9.0 11 123 11.4 4.0 12 124 18.0 7.3 13 111 13.9 4.1 14 124 29.3 6.5 15 121 27.4 9.0 16 118 65.7 13.8 17 122 11.7 3.1 18 123 13.0 4.7 19 125 30.7 4.2 20 133 30.9 6.3 21 122 42.5 5.7 22 123 65.6 6.3 23 119 42.3 4.7 24 119 31.9 4.5 25 120 37.8 4.5 26 116 43.2 7.3 27 115 57.3 8.2 28 130 34.5 5.0 29 114 42.6 6.6 30 113 20.0 2.7 Average 32.9 6.2

TABLE 12 FBPC MPC τ_(I) τ_(D) Total Total Subject K_(c) (min) (min⁻¹) Insulin (U) J Insulin (U) 1 −0.00004 100 600 59.3 600 75.3 2 −0.00004 50 600 77.6 350 94.2 3 −0.00006 100 300 44.7 400 57.8 4 −0.00005 90 450 115.5 600 128.1 5 −0.00006 100 500 94.9 250 113.7

As outlined in Tables 11A and 11B, and as illustrated in FIGS. 11A-E, the control run results of FBPC are considerately better than MPC. In FIGS. 11A-E, FBPC has tighter variation of BGC than MPC in all cases, and much better in some cases (e.g., Subject 5). As given in Tables 11A and 11B, FBPC has smaller Stdev and larger % in the 70-180 mg/dL range, except when both give 100%, than MPC in all 30 cases. Relative to the Stdev of FBPC, the % difference of MPC is as high as 67.2%, greater than 40% for 10 of the 30 subjects, and has an average value of nearly 33% higher. In absolute magnitude, the Stdev of MPC is greater than FBPC; as high as 13.8 mg/dL, greater than 7 mg/dL for 10 subjects, and has an average value of 6.2 mg/dL.

Table 12 gives the values of the tuning parameters for only the first five subjects for space consideration. The values for FBPC are surprising very similar, which may assist greatly in obtaining tuning parameters for FBPC. The values of the one tuning parameter for MPC, J, range from 250 to 600. These values, and the spread between them, are considerably large. The large values are not surprising when considering the large values of the lags shown in Table 10; τ_(i) values for carbohydrates and IFR have averages of 717.3 min. and 920.7 min., respectively.

Table 12 shows the amount of total bolus insulin infused during the control runs of five days. In each case, this amount is significantly higher for MPC. The bolus profiles over time for these five subjects can be seen in FIGS. 12A-E. As shown, for FPBC, the pulses of bolus insulin are more frequent and much smaller than MPC in all cases. For Subject 5, FIG. 12F shows one meal in greater detail. As illustrated, bolus injection for FBPC covers a wider timeframe, is fairly uniform, and its values are considerably lower than MPC, which is also fairly uniform, but less frequent.

Therefore, methods and systems have been described relating to insulin delivery. It should be understood that the present invention contemplates numerous variations, options, and alternatives. For example, where sensed information is provided as input, the sensed information may be supplied from any number of sensors including sensors associated with wearable devices, mobile devices, or other computing devices. Where information is input by users, it is to be understood that the information may be input using wearable devices, mobile devices, or other computing devices including the same device used for insulin delivery. It is further to be understood that the model may be stored on any number of computing devices including a mobile device, or the same device used for insulin delivery.

The above specification provides a description of the manufacture and use of the disclosed compositions and methods. Since many embodiments can be made without departing from the spirit and scope of the invention, the invention resides in the claims. 

What is claimed is:
 1. A method comprising steps of: providing an insulin delivery system; determining the amount of insulin to administer using a feedback predictive control, wherein the feedback predictive control is stored on a machine readable non-transitory media; wherein the feedback predictive control comprises a first model and a second model; wherein the first model is a predictive model; and wherein the second model is a noise model, which models unmeasured disturbances and bias; providing predictive inputs and measured inputs to the first model; parameterizing the first model using the inputs to provide an output, wherein the output is provided to the second model and to a PID controller, wherein the PID controller controls an insulin flow rate; providing feedback to the PID controller from the second model; administering insulin using the insulin delivery system based on the output and the feedback, wherein the administration of insulin is adjusted to nullify the effect of any disturbances on blood glucose concentration based on a predicted blood glucose value.
 2. The method of claim 1 wherein one or more of said measured inputs are provided by an automatic monitoring system.
 3. The method of claim 1 wherein one or more of said predictive inputs are provided manually by a user.
 4. The method of claim 1 wherein the automatic insulin system is a wearable device, an artificial pancreas, or a combination thereof.
 5. The method of claim 2 wherein the automatic monitoring system comprises one or more sensors comprising at least one of the following a soft sensor, a remote sensor, an accelerometer, or a thermistor.
 6. The method of claim 2 wherein the automatic monitoring system monitors at least one of the following variables body position, movement, heat dissipated, skin temperature, near body temperature, galvanic skin response, and sleep, basal insulin, or bolus insulin.
 7. The method of claim 3 wherein the one or more inputs provided manually by the user includes at least one of the following consumed energy, basal insulin, or bolus insulin.
 8. The method of claim 1 wherein the insulin delivery system comprises an apparatus for administering insulin.
 9. The method of claim 1, wherein the first model is: ƒ(V)=η_(t)=ε_(i) =a ₀ +V _(FI,t) +a _(A1) v _(A) _(t) _(,t) + . . . +a _(Ap) v _(A) _(p) _(,t) and wherein the second model is $\begin{matrix} {{\hat{y}}_{t + {k\; \Delta \; t}} = {{\hat{\eta}}_{t + {k\; \Delta \; t}} + {{\hat{\phi}}_{1}\left( {y_{t} - {\hat{\eta}}_{t}} \right)} + {{\hat{\phi}}_{2}\left( {{\hat{y}}_{t - {\Delta \; t}} - {\hat{\eta}}_{t - {\Delta \; t}}} \right)} + \ldots}} \\ {= {{\hat{\eta}}_{t + {k\; \Delta \; t}} + {{\hat{\phi}}_{1}e_{t}} + {{\hat{\phi}}_{2}e_{t - {\Delta \; t}}} + \ldots}} \end{matrix}$ wherein e_(t) is a residual at t of the first model, and subject to: ${{subject}\mspace{14mu} {to}\mspace{14mu} {\sum\limits_{i = 1}^{4}{\hat{\phi}}_{i}}} = 1.$
 10. A system for administering insulin, the system comprising: an insulin delivery system comprising an apparatus for administering insulin; a plurality of sensors; a PID controller comprising a feedback predictive control stored on a machine readable non-transitory media associated with a computing device, the computing device in operative communication with the plurality of sensors; wherein the feedback predictive control comprises a first model and a second model; wherein the first model is a predictive model; and wherein the second model is a noise model, which models unmeasured disturbances and bias; wherein the machine readable non-transitory media is capable of receiving predictive inputs and measured inputs and wherein the first model parameterizes said predictive inputs and/or measured inputs to provide an output, wherein the output is provided to the second model and to a PID controller, wherein the PID controller controls an insulin flow rate, wherein the second model provides feedback to the PID controller, wherein the PID controller controls an insulin flow rate.
 11. The system of claim 10 further comprising a monitoring system, wherein at least one of said measured inputs are provided by the monitoring system.
 12. The system of claim 10 wherein at least one of said predictive inputs are provided manually by a user.
 13. The system of claim 11 wherein the monitoring system monitors at least one of the following variables body position, movement, heat dissipated, skin temperature, near body temperature, galvanic skin response, and sleep, basal insulin, or bolus insulin.
 14. The system of claim 13 wherein the one or more inputs provided manually by the user includes at least one of the following consumed energy, basal insulin, or bolus insulin.
 15. The system of claim 10 wherein the plurality of sensors comprise at least one of the following a soft sensor, a remote sensor, an accelerometer, or a thermistor.
 16. The system of claim 10 wherein the insulin delivery system is a wearable device, an artificial pancreas, or a combination thereof.
 17. The system of claim 10 wherein the first model is: ƒ(V)=η_(t)=ε_(i) =a ₀ +V _(FI,t) +a _(A1) v _(A) _(t) _(,t) + . . . +a _(Ap) v _(A) _(p) _(,t) and wherein the second model is $\begin{matrix} {{\hat{y}}_{t + {k\; \Delta \; t}} = {{\hat{\eta}}_{t + {k\; \Delta \; t}} + {{\hat{\phi}}_{1}\left( {y_{t} - {\hat{\eta}}_{t}} \right)} + {{\hat{\phi}}_{2}\left( {{\hat{y}}_{t - {\Delta \; t}} - {\hat{\eta}}_{t - {\Delta \; t}}} \right)} + \ldots}} \\ {= {{\hat{\eta}}_{t + {k\; \Delta \; t}} + {{\hat{\phi}}_{1}e_{t}} + {{\hat{\phi}}_{2}e_{t - {\Delta \; t}}} + \ldots}} \end{matrix}$ wherein e_(t) is a residual at t of the first model, and subject to: ${{subject}\mspace{14mu} {to}\mspace{14mu} {\sum\limits_{i = 1}^{4}{\hat{\phi}}_{i}}} = 1.$
 18. A system for administering insulin, the system comprising: an insulin delivery system comprising an apparatus for administering insulin; a plurality of sensors comprising at least two of the following a soft sensor, a remote sensor, an accelerometer, or a thermistor; a monitoring system, wherein the monitoring system monitors at least two of the following variables: body position, movement, heat dissipated, skin temperature, near body temperature, galvanic skin response, and sleep, basal insulin, or bolus insulin; wherein the monitoring system is configured to receive one or more inputs provided manually by the user including at least one of the following consumed energy, basal insulin, or bolus insulin; and wherein the plurality of sensors are in communication with the monitoring system; a model stored on a machine readable non-transitory media associated with a computing device, the computing device in operative communication with the plurality of sensors; a PID controller comprising a feedback predictive control stored on a machine readable non-transitory media associated with a computing device, the computing device in operative communication with the plurality of sensors; wherein the feedback predictive control comprises a first model and a second model; wherein the first model is a predictive model; and wherein the second model is a noise model, which models unmeasured disturbances and bias; wherein the machine readable non-transitory media is capable of receiving predictive inputs and measured inputs and wherein the first model parameterizes said predictive inputs and/or measured inputs to provide an output, wherein the output is provided to the second model and to a PID controller, wherein the PID controller controls an insulin flow rate, wherein the second model provides feedback to the PID controller, wherein the PID controller controls an insulin flow rate.
 19. The system of claim 18 wherein the insulin delivery system is a wearable device, an artificial pancreas, or a combination thereof.
 20. The system of claim 19 wherein the first model is: ƒ(V)=η_(t)=ε_(i) =a ₀ +V _(FI,t) +a _(A1) v _(A) _(t) _(,t) + . . . +a _(Ap) v _(A) _(p) _(,t) and wherein the second model is $\begin{matrix} {{\hat{y}}_{t + {k\; \Delta \; t}} = {{\hat{\eta}}_{t + {k\; \Delta \; t}} + {{\hat{\phi}}_{1}\left( {y_{t} - {\hat{\eta}}_{t}} \right)} + {{\hat{\phi}}_{2}\left( {{\hat{y}}_{t - {\Delta \; t}} - {\hat{\eta}}_{t - {\Delta \; t}}} \right)} + \ldots}} \\ {= {{\hat{\eta}}_{t + {k\; \Delta \; t}} + {{\hat{\phi}}_{1}e_{t}} + {{\hat{\phi}}_{2}e_{t - {\Delta \; t}}} + \ldots}} \end{matrix}$ wherein e_(t) is a residual at t of the first model, and subject to: ${{subject}\mspace{14mu} {to}\mspace{14mu} {\sum\limits_{i = 1}^{4}{\hat{\phi}}_{i}}} = 1.$ 